// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_SPARSE_MARKET_IO_H
#define EIGEN_SPARSE_MARKET_IO_H

#include <iostream>
#include <vector>

namespace Eigen {

namespace internal {
template<typename Scalar, typename StorageIndex>
inline void
GetMarketLine(const char* line, StorageIndex& i, StorageIndex& j, Scalar& value)
{
	std::stringstream sline(line);
	sline >> i >> j >> value;
}

template<>
inline void
GetMarketLine(const char* line, int& i, int& j, float& value)
{
	std::sscanf(line, "%d %d %g", &i, &j, &value);
}

template<>
inline void
GetMarketLine(const char* line, int& i, int& j, double& value)
{
	std::sscanf(line, "%d %d %lg", &i, &j, &value);
}

template<>
inline void
GetMarketLine(const char* line, int& i, int& j, std::complex<float>& value)
{
	std::sscanf(line, "%d %d %g %g", &i, &j, &numext::real_ref(value), &numext::imag_ref(value));
}

template<>
inline void
GetMarketLine(const char* line, int& i, int& j, std::complex<double>& value)
{
	std::sscanf(line, "%d %d %lg %lg", &i, &j, &numext::real_ref(value), &numext::imag_ref(value));
}

template<typename Scalar, typename StorageIndex>
inline void
GetMarketLine(const char* line, StorageIndex& i, StorageIndex& j, std::complex<Scalar>& value)
{
	std::stringstream sline(line);
	Scalar valR, valI;
	sline >> i >> j >> valR >> valI;
	value = std::complex<Scalar>(valR, valI);
}

template<typename RealScalar>
inline void
GetVectorElt(const std::string& line, RealScalar& val)
{
	std::istringstream newline(line);
	newline >> val;
}

template<typename RealScalar>
inline void
GetVectorElt(const std::string& line, std::complex<RealScalar>& val)
{
	RealScalar valR, valI;
	std::istringstream newline(line);
	newline >> valR >> valI;
	val = std::complex<RealScalar>(valR, valI);
}

template<typename Scalar>
inline void
putMarketHeader(std::string& header, int sym)
{
	header = "%%MatrixMarket matrix coordinate ";
	if (internal::is_same<Scalar, std::complex<float>>::value ||
		internal::is_same<Scalar, std::complex<double>>::value) {
		header += " complex";
		if (sym == Symmetric)
			header += " symmetric";
		else if (sym == SelfAdjoint)
			header += " Hermitian";
		else
			header += " general";
	} else {
		header += " real";
		if (sym == Symmetric)
			header += " symmetric";
		else
			header += " general";
	}
}

template<typename Scalar, typename StorageIndex>
inline void
PutMatrixElt(Scalar value, StorageIndex row, StorageIndex col, std::ofstream& out)
{
	out << row << " " << col << " " << value << "\n";
}
template<typename Scalar, typename StorageIndex>
inline void
PutMatrixElt(std::complex<Scalar> value, StorageIndex row, StorageIndex col, std::ofstream& out)
{
	out << row << " " << col << " " << value.real() << " " << value.imag() << "\n";
}

template<typename Scalar>
inline void
putVectorElt(Scalar value, std::ofstream& out)
{
	out << value << "\n";
}
template<typename Scalar>
inline void
putVectorElt(std::complex<Scalar> value, std::ofstream& out)
{
	out << value.real() << " " << value.imag() << "\n";
}

} // end namespace internal

inline bool
getMarketHeader(const std::string& filename, int& sym, bool& iscomplex, bool& isvector)
{
	sym = 0;
	iscomplex = false;
	isvector = false;
	std::ifstream in(filename.c_str(), std::ios::in);
	if (!in)
		return false;

	std::string line;
	// The matrix header is always the first line in the file
	std::getline(in, line);
	eigen_assert(in.good());

	std::stringstream fmtline(line);
	std::string substr[5];
	fmtline >> substr[0] >> substr[1] >> substr[2] >> substr[3] >> substr[4];
	if (substr[2].compare("array") == 0)
		isvector = true;
	if (substr[3].compare("complex") == 0)
		iscomplex = true;
	if (substr[4].compare("symmetric") == 0)
		sym = Symmetric;
	else if (substr[4].compare("Hermitian") == 0)
		sym = SelfAdjoint;

	return true;
}

template<typename SparseMatrixType>
bool
loadMarket(SparseMatrixType& mat, const std::string& filename)
{
	typedef typename SparseMatrixType::Scalar Scalar;
	typedef typename SparseMatrixType::StorageIndex StorageIndex;
	std::ifstream input(filename.c_str(), std::ios::in);
	if (!input)
		return false;

	char rdbuffer[4096];
	input.rdbuf()->pubsetbuf(rdbuffer, 4096);

	const int maxBuffersize = 2048;
	char buffer[maxBuffersize];

	bool readsizes = false;

	typedef Triplet<Scalar, StorageIndex> T;
	std::vector<T> elements;

	Index M(-1), N(-1), NNZ(-1);
	Index count = 0;
	while (input.getline(buffer, maxBuffersize)) {
		// skip comments
		// NOTE An appropriate test should be done on the header to get the  symmetry
		if (buffer[0] == '%')
			continue;

		if (!readsizes) {
			std::stringstream line(buffer);
			line >> M >> N >> NNZ;
			if (M > 0 && N > 0) {
				readsizes = true;
				mat.resize(M, N);
				mat.reserve(NNZ);
			}
		} else {
			StorageIndex i(-1), j(-1);
			Scalar value;
			internal::GetMarketLine(buffer, i, j, value);

			i--;
			j--;
			if (i >= 0 && j >= 0 && i < M && j < N) {
				++count;
				elements.push_back(T(i, j, value));
			} else
				std::cerr << "Invalid read: " << i << "," << j << "\n";
		}
	}

	mat.setFromTriplets(elements.begin(), elements.end());
	if (count != NNZ)
		std::cerr << count << "!=" << NNZ << "\n";

	input.close();
	return true;
}

template<typename VectorType>
bool
loadMarketVector(VectorType& vec, const std::string& filename)
{
	typedef typename VectorType::Scalar Scalar;
	std::ifstream in(filename.c_str(), std::ios::in);
	if (!in)
		return false;

	std::string line;
	int n(0), col(0);
	do { // Skip comments
		std::getline(in, line);
		eigen_assert(in.good());
	} while (line[0] == '%');
	std::istringstream newline(line);
	newline >> n >> col;
	eigen_assert(n > 0 && col > 0);
	vec.resize(n);
	int i = 0;
	Scalar value;
	while (std::getline(in, line) && (i < n)) {
		internal::GetVectorElt(line, value);
		vec(i++) = value;
	}
	in.close();
	if (i != n) {
		std::cerr << "Unable to read all elements from file " << filename << "\n";
		return false;
	}
	return true;
}

template<typename SparseMatrixType>
bool
saveMarket(const SparseMatrixType& mat, const std::string& filename, int sym = 0)
{
	typedef typename SparseMatrixType::Scalar Scalar;
	typedef typename SparseMatrixType::RealScalar RealScalar;
	std::ofstream out(filename.c_str(), std::ios::out);
	if (!out)
		return false;

	out.flags(std::ios_base::scientific);
	out.precision(std::numeric_limits<RealScalar>::digits10 + 2);
	std::string header;
	internal::putMarketHeader<Scalar>(header, sym);
	out << header << std::endl;
	out << mat.rows() << " " << mat.cols() << " " << mat.nonZeros() << "\n";
	int count = 0;
	for (int j = 0; j < mat.outerSize(); ++j)
		for (typename SparseMatrixType::InnerIterator it(mat, j); it; ++it) {
			++count;
			internal::PutMatrixElt(it.value(), it.row() + 1, it.col() + 1, out);
		}
	out.close();
	return true;
}

template<typename VectorType>
bool
saveMarketVector(const VectorType& vec, const std::string& filename)
{
	typedef typename VectorType::Scalar Scalar;
	typedef typename VectorType::RealScalar RealScalar;
	std::ofstream out(filename.c_str(), std::ios::out);
	if (!out)
		return false;

	out.flags(std::ios_base::scientific);
	out.precision(std::numeric_limits<RealScalar>::digits10 + 2);
	if (internal::is_same<Scalar, std::complex<float>>::value || internal::is_same<Scalar, std::complex<double>>::value)
		out << "%%MatrixMarket matrix array complex general\n";
	else
		out << "%%MatrixMarket matrix array real general\n";
	out << vec.size() << " " << 1 << "\n";
	for (int i = 0; i < vec.size(); i++) {
		internal::putVectorElt(vec(i), out);
	}
	out.close();
	return true;
}

} // end namespace Eigen

#endif // EIGEN_SPARSE_MARKET_IO_H
